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Abstract We investigate the use of Malliavin calculus in order to calculate the Greeks of 
multidimensional complex path-dependent options by simulation. For this purpose, we ex- 
tend the formulas employed by Montero and Kohatsu-Higa to the multidimensional case. 
The multidimensional setting shows the convenience of the Malliavin Calculus approach 
over different techniques that have been previously proposed. Indeed, these techniques may 
be computationally expensive and do not provide flexibility for variance reduction. In con- 
trast, the Malliavin approach exhibits a higher flexibility by providing a class of functions 
that return the same expected value (the Greek) with different accuracies. This versatility for 
variance reduction is not possible without the use of the generalized integral by part formula 
of Malliavin Calculus. In the multidimensional context, we find convenient formulas that 
permit to improve the localization technique, introduced in Fournie et al and reduce both 
the computational cost and the variance. Moreover, we show that the parameters employed 
for variance reduction can be obtained on the flight in the simulation. We illustrate the ef- 
ficiency of the proposed procedures, coupled with the enhanced version of Quasi-Monte 
Carlo simulations as discussed in Sabino, for the numerical estimation of the Deltas of call, 
digital Asian-style and Exotic basket options with a fixed and a floating strike price in a 
multidimensional Black-Scholes market. 

Key Words: Greeks, Risk-Management, Quasi-Monte Carlo Methods, Malliavin Calculus. 



1 Introduction and Motivation 

Risk-sensitivities, also called Greeks, are fundamental quantities for the risk-management. 
Greeks measure the sensitivities of a portfolio of financial instruments with respect to the 
parameters of the underlying model. Mathematically speaking, a greek is the derivative of 
a financial quantity with respect to (w.r.t.) any of the parameters of the problem. As these 
quantities measure risk, it is important to calculate them quickly and with a small order of 
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error. In general, the computational effort required for an accurate calculation of sensitivities 
is often substantially greater than that required for price estimation. 

The problem of greeks calculation can be casted as follows. Suppose that the financial 
quantity of interest is described by E[\f/{X{a))Y\ (i.e., the price of a derivative contract), 
where t|/ : M — > R is a measurable function and X and Y are two random variables (r.v.s). 
The greek, that we denote 6, is the derivative w.r.t. the parameter a : 



e(a) = ^E[v/(X(a))F]=] 



4-¥{X{a))Y 
da 



The most common of the Greeks are notably. Delta, Gamma, Vega, Theta, Rho. These quan- 
tities are relatively simple to calculate for plain vanilla contracts in the Black-Scholes (BS) 
market. However, their evaluation is a complex and demanding task for exotic derivative 
contracts such as Asian-style basket options where no closed-formula is known. 

The simplest and crudest approach is to employ the Monte Carlo (MC) estimation of 
¥i[\if{X{a))Y] for two or more values of a and then use finite-difference approximations. 
However, this approach can be computationally intensive and can produce large biases and 
large variances in particular if y/^ = 11^, where A is a measurable set. A variant is the kernel 
method (see Montero and Kohatsu-Higa IIOII ) which generalizes finite-difference methods 
using ideas taken from the kernel density estimation. 

Several alternatives have been proposed without finite-difference approximation. Path- 
wise methods (see Glasserman (4)) treat the parameter of differentiation a as a parameter 
of the evolution of the underlying model and differentiate this evolution. However, this ap- 
proach is not always applicable, notably when )|/ is not smooth (for instance )|/ = U^). At 
the other extreme, the likelihood method ratio (see Glasserman |4|) puts the parameter in 
the measure describing the underlying model and differentiates this measure. Even if the 
likelihood method ratio is applicable to non-smooth functions it may provide high- variance 
estimators. Indeed, compared to the pathwise method (when applicable), it displays a higher 
variance. Summarizing, these two alternatives involve two main ideas: differentiating the 
evolution or differentiating the measure, respectively. 

In this paper we investigate the use of Malliavin Calculus in order to employ (Quasi)- 
Monte Carlo (QMC) simulations for the evaluation of the sensitivities of complex multi- 
dimensional path-dependent options. The multidimensional setting shows the very conve- 
nience of the Malliavin Calculus approach over the different techniques that have been pro- 
posed. Indeed, Malliavin Calculus allows to calculate sensitivities as expected values whose 
estimation is a natural application of MC methods. Formally: 



e(a) = E 



^W{Xia))Y 



-■¥.[^{X{a))H] 



where // is a r.v. depending on X and Y . 

In the context of multidimensional options, we extend the formulas employed by Mon- 
tero and Kohatsu-Higa lllll to the multidimensional case. This approach gives a certain 
flexibility and provides a class of functions (different r.v.s H) returning the same expected 
value (the sensitivity) but with different accuracies. 

Indeed, the previously mentioned alternative techniques may be computationally expen- 
sive in the multidimensional case and do not provide flexibility for variance reduction. This 
versatility for variance reduction is not possible without the use of the generalized inte- 
gral by part formula of Malliavin Calculus. Advanced techniques such as the kernel density 
estimation or more recent approaches such as the Vibrato Monte Carlo in Gilles |3| are 
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difficult to employ and computationally demanding in multi-dimensions. In order to avoid 
to use the Malliavin technique, Chen and Glasserman 1 1] have illustrated a procedure that 
produces "Malliavin Greeks" without Malliavin Calculus. However, since this procedure in- 
volves both pathwise and likelihood ratio methods, the estimators of the formulas for the 
sensitivities in Chen and Glasserman [ 1] have a high variance. 

For these purposes we find convenient representations of H that permit to enhance the 
localization technique introduced in Foumie et al. |2| and reduce both the computational 
cost and the variance. Moreover, we show that the parameters employed for the variance 
reduction can be obtained on the flight in the simulation by adaptive techniques. We illus- 
trate the efficiency of the proposed procedures, coupled with the enhanced version of QMC 
simulations discussed in Sabino |16|, for the numerical estimation of the Deltas of call, 
digital Asian-style and Exotic basket options with a fixed and a floating strike price in a 
multidimensional BS market. 

The paper is organized as follows. Section |2] is a short introduction on Malliavin Cal- 
culus, Section [3] derives the formulas employed for the computation of the Deltas of call 
Asian basket options with floating and fixed strike, Asian digital options and exotic options. 
Section|4]illustrates the enhanced QMC approach that we adopt and describes in details how 
to get the localization parameters with adaptive (Q)MC techniques; Section|5]discusses the 
numerical experiments of the study and finally Section |6] summarizes the most important 
results and concludes the paper. 

2 Malliavin Calculus: Basic Results and Notation 

The aim of this section is to briefly introduce the basic results from Malliavin Calculus and 
to fix the notation we adopt in the rest of the paper. For more information on this subject, 
we refer the reader to the book by Nualart 11131 . 

Consider the probability space (12, ^,P) where we define the M-dimensional Brown- 

ian motion W(r) = (Wi(r), . . . , Wm(0), ? e [0,r] and given = . . . ,4"', divide the 

interval [0, T] into n subintervals 4 = , /^"'j ),fc = 0,...,n— 1. The superscripts indicates 
the fineness of the subdivision of [0,7^]. Now denote the vector 4^"'' = ^Zi]"' , . . . , Zi^']||, j 

where 4<;:.'=Ty,(4:\)-^-(4"')- 

Now consider a smooth function with polynomial growth ^ : R"^*' — )• R, 6 of 

the form (j) = (j) (4q"' , . . . , ) . Finally we consider the following space: 

^, = { (4 . . . , Ai"\ y,^e'^;}c^\Q), (1) 
where (.-$1,)„>i form an increasing sequence in Jif^ {Q). 

Definition 1 The union 5^ = (lj«>i •S^n) C iQ.) is called simple functional space and 
its elements are called simple functionals. 

We now can define the Malliavin derivative operator. 

Definition 2 Let (j) e ,y, then there exists « 6 N* such that (j) = (j){A^"^). The Malliavin 
derivative operator D = (Z)' , . . . ,D^^ of(j) at a point s 6 If. is defined as 

^"'<?=E^(^o"\-->4i"\)U^,„,(.). (2) 
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Let us precise the notation. We have xj^ = (xj;.i , . . . ,x,t jvf), so xj; corresponds to the increment 

in) in) 

vector Aj^ and the m-th component x^^m corresponds to the m-th component Aj^ ^. 
Moreover we give the following definition. 
Definition 3 Introduce on ,y the norm 

\\F\\\.2 = \\F\\.^^-(n) + l|D^'||^2([o,7-|xfi)' 

the ief D''^ is the closure of 5^ with respect (w.rt.) \\ ■ |li_2- 

Finally we define tlie Skorohod integral 5^^. 

Definition 4 The adjoint o/D in {Q x [0, T]) is the operator: 



5'* : u = (mi , . . . , M„) 6 dom (5'*) ^ 5^* ( 



which by definition satisfies for 6 D''^ and u e dom (5'^*) 



M 

EE 

m— 1 



Df {(j)u,„{s))ds 



m=l 



:E[(^5^*(u) 



(3) 



(4) 



Equation (O is known as duality relation. 

It can be shown (see Nualart fTTl) that if u(f) is an Ito process, the Skorohod integral 
coincides with the Ito integral of u and DjU = if i > f . 

We now list some identities and useful results that will be employed in the rest of this 
paper. Proofs can be found in Nualart Ill3i . 

1. VFi,...,Fj eD''2 we have (j) (Fu . . . , Fj) eD''^ and Vm = 1,...,M 



D'^<t> = Z^^{Fu-.;F,)D:F,. 



(5) 



For example, let a e R'" we have 

' M 



D";expl Y^atWiit) =a,„exp J^a,W,(f) U[o,,](j 



2. Let (j)(t) be an adapted process we have: 
rT 



!=1 



D' 







(l){t)dW"'{t) = (l){s) + / Df(l){t)dW'"{t) 



and 



(l){t)dt = 1^ D";^{t)dt 
3. If e B^'^ u e dom (5^''), and (^u e dom {8^^) then 



M /.J 



5S''((?u)=(?5^''(u)- / «,„(.)D7(?)J^. 



(6) 



(7) 



(8) 



(9) 
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3 Multidimensional Malliavin Sensitivities 

Consider for simplicity a complete market whose risky assets, S,-, / = 1, . . . ,M, are driven by 
the following dynamics (in the risk-neutral measure): 

dSi(t) = rSi{t)d[ + Si{t)ai{t)dBi{t) i=l,...,M, (10) 
Si{0) = Xi, 

where r is the constant risk-free rate, a{t) = {oi{t), . . ., cTm(O) '^he vector of the volatilities 
process and B(f) is the vector of the M-dimensional Brownian motion in the risk-neutral 
measure with dBi{t)dB„i{t) = Pi,„{t)dt; p is the correlation matrix among the Brownian 
motions (it can be stochastic). The existence of the vector process a{t) is guaranteed by 
theorem 9.2.1 in Shreve [181. Applying the risk-neutral pricing formula (see Shreve [18.1 ), 
the calculation of the price at time / of any European derivative contract with maturity date 
T boils down to the evaluation of an (discounted) expectation: 

a{t)=exp{-r{T-t))E[xif\^,], (11) 

the expectation is under the risk-neutral probability measure and xj/isa generic .^j- -measurable 
variable that determines the payoff of the contract. 

In order to apply Malliavin Calculus we need to write the above dynamics in terms of 
uncorrelated Brownian motions: 

M 

dSi{t) = rSi{t)dt + Si{t)ai{t) ai,n{t)dW„,{t) i=l,...,M 



m— 1 



5,(0) = Xi, 



where Y.m^^ a,>„(r)afa„(r) = Pik{t),a.s. and we have defined CT,„,(f) = Oi{t)Y.m^^ aim{t),a.s.. 

Hereafter we denote 5*^' and 8^ the Kronecker delta and the Dirac delta, respectively. 
Naturally at time T we have )|/ = a{T), a.s.. 

The following proposition generalizes the formula in Montero Kohatsu-Higa lllll to the 
multidimensional case. 

Proposition 1 Assuming the dynamics ( I72t let m(T) be a -measurable r.v. ( it can depend 
on the entire trajectory) and consider \ff = \ff{m{T)). Denote Gk the partial derivative 

dm(T) 

G, = ^^, k=l,...,M, (12) 

OXk 



Suppose that iff 6 D''^, the k-th delta (the k-th component of the gradient) is 



M 

¥ E ^miGkUm) 



(13) 



m— 1 

where u = (mi, . . .,um) £ dom(5^''), z = (zi, . . .,Zm) G dom(5^''), G^u e dom(5^'') and 

^ / N ^ z,n{s) 

LtJo Zhis)D^m{T)ds 

M „T 

£ / Zh{s)D'MT)ds ^ 0, a.s. 
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The derivative may have no mathematical sense indeed, the aim of the proposition is to 
overtake the problem with the formalism of distributions and Malliavin Calculus. 



Proof Compute 



D'^\ff=\f/D'lm{T) h=\,...,M. 



(14) 



Suppose z e dom(5^'^) and multiply the above equation by z/,(r) and by Gk; then sum for all 
h = 1 , . . . , M and integrate: 



£/ G,Zh{s)D''MT)ds = GkZh{sW{T)D'lm{T)ds. (15) 



M nj 



y/'Gk does not depend on s and due to the definition of u we can write 



M pT 



;))— 1 



iff'Gk = / u,„{s)GkDfiff{T))ds. k=\,...,M 



(16) 



Finally compute the expected value of both sides of ( 116b 
E [y/Gk] = E 

By duality 



M nT 

/ u„{s)GkDfxifds 



Ak=E^\j/5^^{Gkyi)\ k=\.....,M, 
and this concludes the proof. 



(17) 



(18) 



3.1 Greeks in the Multidimensional Black-Scholes Market 

In this section we apply Proposition [T] to the case of a multidimensional Black-Scholes 
market where the volatilities vector process in Equation (112b is not stochastic (for simplicity 
we consider constant volatilities and correlations). The main advantage of the Malliavin 
approach over different techniques, for example the methods in Gilles |3 1 and the Chen and 
Glasserman H I, is that Proposition [T] allows the possibility of variance and computational 
reduction due to the flexibility in choosing either the process u, or better z. The methods 
illustrated in Gilles |3 | and Chen and Glasserman |[T] are difficult to employ if we assume a 
multidimensional dynamics and they do not allow versatility for variance reduction. 

We consider the case Zh = (^kSj^', h,k = I,. . . ,M, ak = l,Vfc. Namely, in order to com- 
pute the k-th delta we consider only the k-th term of the Skorohod integral reducing the 
computational cost. In particular, this choice is motivated by the fact that we can enhance 
the localization technique introduced by Foume et al. [2|. With this setting we need to con- 
trol only 8^^{-) and then only the A;-th component of W(f). This enhancement is not possible 
with other approaches that furnish only a fixed representation of the components of the mul- 
tidimensional deltas. 

Under the above assumptions for the vector process z, we explicitly derive the multidi- 
mensional deltas for the following exotic options in the BS market: 
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1. Discretely monitored Asian basket options with fixed strike. Assume fi < t2 ■ ■ ■ < = T , 
where T is the maturity of the contract and the payoff function 

/' M N \ + 

¥=iLL'''uSiitj)-Kj , (19) 

where K is the strike price and Wij = 1 . In this case we have 

1 

Gk = — "^kjSkitj) 



and 



M N 

We then calculate the following quantities 

Lk= / D'',m{T)ds = E E ^ijSi{tj)tj^ik^ 



Ak = 



,=1 j=i 

T N N 

D'lGkds=Y^WjkSk{tj)tjakk{s) = E WjkSk{tj)tjak 

r N 
D%ds=Y^WijSi{tj)t]afk, 



and hence 



4t = E 



Sk 



, /t=l,...,M. 



(20) 



Due to the equation ^ we can write the the Skorohod integral above for k= 1, . . . ,M 
as: 



Lk J Lk Li 



DlGkds-Gk DlLkds 



Lk 



Wk{T) + - 



Bk\ Ak 



(21) 

With another choice of z, for instance Zh = cCh, would depend linearly on the whole 
M-dimensional Brownian motion, making the localization technique less efficient. 

2. Discretely monitored Asian basket options with floating strike K{T) - 



M 



-. For 

simplicity we assume Wij = j^^li,]. The calculation is similar to the previous payoff 
function, indeed we can write \j/ = n{T)^ where 

n{T)=m{T)-K{T). 

In analogy, we have 

dn(T) Sk(T) 
rk = = Gk — = Gk — Ik, 



Mk 



Mxk 



D,n{T)ds = Lk- D,K{T) = Lk 



LfLiSi{T)Taik 



M 



■ Lk — Uk, 



r D'lFkds^Ak - f D%ds = Ak ■ 
Jo Jo 



Sk{T)Tak 



Mxk 



--Ak-Vk, 



DlMkds= D^Lkds- DlUkds = B, 



ltiSi(T)T^atk 



2„2 



M 



- Bk —Pk 
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with the quantities Tk,Uk;Vk,Pk,^k automatically defined by the above equations. Then 



and 



4t = E 



El 

Mk 



, k=l,...,M. 



MkJ Mk 



Wk{T) + -- 



Bk-Pk\ Ak-Vk 



Mk 



Mk 



3. Digital Asian basket options with fixed strike. 

¥= ^m(T)>K- 



(22) 



(23) 



(24) 



This type of payoff function fulfills the hypotheses of Proposition[T]and we might adopt 
equation ilOi . However, due to the properties of the Dirac delta 5^ and Proposition [T] 
we can write 



Ak = e"''^E [5^ {m{T))Gk] = e-''^E 



i{T)-K 



Gk 



M 



¥ E S^^{Gk(l)u^) 



where we assume that (j), (j)' are square integrable, (j){0) = 1, (j)Gk is Skorohod integrable 
Vfc = 1, . . . ,M and h>0. The aim of this setting is to reduce the variance of the MC esti- 
mator of All by tuning the localization function (j) around the strike K with a convenient 
choice of the parameter h (see Kohatsu-Higa and Patterson |7|). 
Under this assumption the Skorohod integral in equation (I13l l becomes: 



GkU 



m{T)-K 



M i-T 



D7 U 



<T)-K 

then the Skorohod integral 5^^ 

'm(T)-K 



<T)-K 
h 

m(T)-K 



Giu) is 



Gkd^\vi) 



m(T)-K \ 
h 



Lm=i lo ^mi'W'Gkds Gk ,fm{T)-K 



LTn=Jo '*m{s)Dfm{T)ds 



Finally, with our choice for the simple process u the last equation becomes 

^m{T)-K' 



i{T)-K 



Lk 



^_(h.,(m{T)-K 
h ^ 



(25) 



(26) 



where 5^''(u) depends on the terms that we have found in the case of Call Asian basket 
options. 

It is worthwile to say that the same localization procedure and the Malliavin approach 
adopted for digital options can be employed for the computation the Gamma (second 
order derivative) for Call Asian basket options. 
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3.2 Greeks for Exotic Options 

In Proposition [T] we have supposed that the payoff function \j/ depends on m{T) only. With 
the notation adopted in the BS setting, suppose for instance, that \j/ = max {m{T) — K, K{T) — K,0), 
where is a fixed price, now we cannot rely on Proposition [T] to derive the expression of 
the sensitivities of such an exotic option. Here \j/ depends separately on two random vari- 
ables m{T) and K{T). In the following we extend Proposition [T] in order to allow such a 
dependence. 



Proposition 2 Assuming the dynamics | |72I ) suppose \j/ = \j/{X,Y). For simplicity we set 
'■^ and Tk = 4^. 



r = Q, denote = ^ and T^- = J^. Let u and p be two simple processes belonging to 



dom(5 ). Define the following -measurable r.v.s: 

M „T 

«i=E/ ^m{s)DfXds, a2 = '£'^,^J^u„,{s)D'^'Yds (27) 

m=l "'0 

^i=E/ Pm{^)DTX, b2 = L^,^i,kJp,n{s)DfYds (28) 

M „T 

Oi = E / nuUs)D'^xifds, 02 = Ijf=, JcJ G,p„{s)Dfxifds (29) 



m— 1 



b2-^ si- 
aib2—a2b\ ' 



^1= -,- : U2 = j^t^. (30) 



Finally, suppose that aib2 — 02^1 7^ 0,a.s. and t/iTiu — (/2GiP is Skorohod integrable, we 
have: 



Ak = =A,=E [T,dxWiX,Y) + GkdYXj/{X,Y)] = E 



V/(X,F) E €(.UlTkU,n-U2G,p,n) 
m— 1 

(31) 

where dx and dy denote the partial derivatives with respect to the first and second variable, 
respectively. 



Proof Compute: 

D'l' \i/{X, Y) = dx WD'I'X + dy \lfD"^Y. (32) 
As done in the proof of Proposition [T] multiply for and m„,, sum for all m and integrate: 



y nu„{s)Dfxifds = y nu„{s)D'^Xds+y Ti,u„{s)DfYds. (33) 

m=l-^0 m=l-^0 

Now repeat the procedure above considering and p/, (s) , we have 



M pT M pT M pT 

E / G,p„(i)D»VJi= E / GkP,n{s)iy':Xds+Y. / GkP,n{s)D^Yds. (34) 

m=l "'0 m=l "'0 m=l "'0 

We rewrite Equations ( I33t and (I33t as a linear system 

Oi= aiTkdx\tf + a2TkdY\lf 
02 = biGkdxW + g2GkdYW 
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Our aim is to compute Tiidxy/{X,Y) + GkdY\l/{X,Y) such that we can apply the duality 
relation. After some algebra we get that 



Gk{aib2-a2bi) ' 

aiTkOi-biGkOi 
Tk{a\b2-a2bi) 



and 



then we have 



Tkdx\lf + GkdY\lf'- 



Ol [b2 



n 



Ak = E[OiUi-02U2]=E 
by duality 



M 

E 

m=I 



01^2 —02^1 



{Ui TkUm [s) - U2GkPm (s) ) Wds 



M 



WY,S^HUlTkU,n-U2GkP„ 



(36) 



(37) 



(38) 



and this concludes the proof. 



We can adapt the result of Proposition |2] to the BS market. Again the Malliavin Calculus 
approach is very versatile and permits to reduce the computational burden and the variance 
of the MC by enhancing the localization technique. As done before we consider (s) = 
5^,\/s and p„,{s) = s5^,\/s, in order to fulfill the hypothesis of Proposition |2] 
The formula for the k-th component of the delta is 



Ak = K [xi/{X, Y) (5f {UiTk) - 5f (^t/zG,))] , 
and the two Skorohod integrals are respectively: 

5f {Ui Tk) = Ui TkWk{T) - Ui r D%ds - o'^.Uids, 

Jo Jo 



(39) 



(40) 



5f {U2Gk) = U2Gk r sdW^ - Gk r sD%ds - U2 f sO'^fikds. (41) 
Jo Jo ' Jo 

In the MC estimation we can simulate the first term in the above equation relying on the 
equality: 

sdW\s) = TWk{T)- Wk{s)ds, 
Jo Jo 

where Jq Wk {s)ds is approximated by a sum at the points t\,. . . ,tN = T . 

In our numerical experiments we consider x^f = max(m(r) —K,K{T) — ^^,0) where 
m{T) and K(T) have been defined in Section [STl The terms in Equations ( 14Gb and i l41b 
have been obtained as in SectionlSTl 
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4 Simulation Setting 

In this section we briefly describe the numerical setting that we adopt for the QMC esti- 
mation of the Greeks by the Malliavin approach formulas. We briefly illustrate the QMC 
method and discuss how to conveniently find the parameters of the localization technique 
on the fly by adaptive simulation. 



4.1 The Quasi-Monte Carlo Framework 

Consider / = E [t|/(X)] where X is a J-dimensional random vector and t|/ : R'' — > R, the QMC 

estimator of / is Iqmc = ""'Jg — ' is the number of simulations, as for the standard 
MC. However the points X, are not pseudo-random but are obtained by low-discrepancy se- 
quences. Low-discrepancy sequences do not mimic randomness but display better regularity 
and distribution (see Niederreiter [12] for more on this subject). We do not enter into the de- 
tails of QMC methods and their properties, we just stress the fact that such techniques do not 
rely on the central limit theorem and the error bounds are given by the well known Hlawka- 
Koksma inequality. Some randomness is then introduced in order to statistically estimate 
the error of the estimation by the sampled variance; this task is achieved by a technique 
called scrambling (see Owen IllSJ ). The randomized version of QMC is called Randomized 
Quasi-Monte Carlo (RQMC). 

In our numerical estimation we use a randomized version of the Sobol' sequence with 
SoboFs property A, that is one of the most used low-discrepancy sequences (it is also a 
digital net). 

Finally, in order to improve the efficiency of RQMC and reduce the effect of the so- 
called curse of dimensionality, we employ the Linear Transformation (LT) technique intro- 
duced in Imai and Tan (5l in the enhanced version illustrated in Sabino I16III7I . The aim 
of the LT algorithm is to concentrate the variance of t// into the components with higher 
variability so that we may profit from the higher regularity of low-discrepancy points and 
then reduce the nominal dimension of 

We briefly describe the LT algorithm. Consider a d dimensional normal random vector 
T ~ ,yV[ix;Z), a vector w = {w\ ,.. ., w^) 6 M.^ and let /(T) = T,i=i ^iTi be a linear com- 

bination of T. Let C be such that E = CC and assume e ~ ,yV{0,Id) with T = Ce. The 
LT approach considers C as C = C^^ = C*^^A, with C'^^ the Cholesky decomposition of E. 
Then, in the linear case, we can define: 



where = C"]^ • w = A.k -B, k= \ ...,d and B = (C'^^)^w while Ck and A.k are the A:-th 
columns of the matrix C and A, respectively. In the linear case, setting 



d 



/(e) :=/(C™Ae) = £ a^e^ + M-w, 



(42) 




(43) 



with arbitrary remaining columns with the only constrain that AA^ 
ing expression: 

/(e)=/i-w±||B||ei. 



= Id, leads to the follow- 
(44) 
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This is equivalent to reduce the effective dimension in the truncation sense to 1 and this 
means to maximize the variance of the first component Ei . 

In a non-linear framework, we can use the LT construction, which relies on the first 
order Taylor expansion of : 

/(e)«/(e) + X:^^4e,. (45) 

The approximated function is linear in the standard normal random vector Ae ~ ■yV{0,Id) 
and we can rely on the considerations above. The first column of the matrix A* is then: 



2 



A.i = arg max — r (46) 

A.16R'' V •'^i / 

Since we have already maximized the variance contribution for ^^^^^ , we might con- 
sider the expansion of g about d — I different points in order to improve the method using 
adequate columns. More precisely Imai and Tan Q propose to maximize: 

A.k = arg max — (47) 

subject to ||A.k*|| = 1 and A.j* ■ A.k* = OJ = I, . . . ,k- l,k < d. 

Although equation ( 143 b provides an easy solution at each step, the correct procedure 
requires that the column vector A.k* is orthogonal to ail the previous (and future) columns. 
Imai and Tan |5| propose to choose e = ei = E[e] = 0, £2 = (1,0, . . . ,0), . . = (1, 1, 1, . . . ,0, . 
where the k-th point has k — I leading ones. We refer to Sabino 11611171 for the details of a 
fast and convenient implementation of this algorithm. 



4.2 Enhancing the Localization Technique 

The aim of the localization technique introduced in Foumie et al. fZ) is to reduce the variance 
of the MC estimator for the sensitivities by localizing the integration by part formula around 
the singularity. In the following, for simplicity, we illustrate the localization technique in the 
case of vanilla call options. 

Foumie et al. I2j found that a (possible) expression for the delta of a call option is: 



A = -r-Ele'^ (S(T)-K)+] =E 
ax 



'iS{T)-Ky 



W{T) 
xTa 



(48) 



When the one-dimensional Brownian motion W{T) is large, the term {S{T) — K)'^ W {T) 
becomes even larger and has a high variance. The idea is to introduce a localization function 
around the singularity at K. 



For 5 > 0, set 



0, 



2S 
1 



for y<K-5, 
for ye [K-d,K + 5], 
for y>K + 5, 



(49) 
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and G{z) = J^aoHg{y)dy, then consider Fg{z) = (z — K)'^ — Gs{z). Consequendy, we have: 



: e' ' E 



HsiS{T)) 



dS{T) 



W{T) 
xTa 



(50) 



Fg vanishes forz <K — 5 and z> K + 5, thus Fg{S{T))W{T) vanishes when W{T) is large. 

The same analysis, with similar results, is valid for the call-style Asian options and the 
exotic option analyzed in Section [3] Indeed, it suffices to replace 5(7') with the average 
jW/y5,(fj) in the equations above and consider an if, else statement to select the local- 
ization function when the strike price is stochastic or the option is exotic. In addition, in 
the above options formulas, the role of the "weight" term ^^^^^ is played by the Skorohod 
integral. We remark that the formulas that we derived to compute the k-th component of the 
delta display weights that depend only on the Skorohod intergral w.r.t. the k-\h component 
of the multidimensional Brownian motion permitting to better control the variance. If we 
would have chosen to control all the components of the Skorohod integral, taking all non- 
zero components of the simple vector process u, we would have needed to tune different 
M Brownian motions making the localization technique less efficient and computationally 
more expansive. 

The choice of the parameter 5 is of fundamental importance for the result of the local- 
ization technique because it influences the variance of the MC estimator. In the following 
we describe how to employ an on the fly efficient value based on adaptive MC simulations. 
For ease of notation, we consider once more a vanilla call option payoff bearing in mind that 
the same applies to the payoffs under study. In such cases we need to make the substitution 
illustrated above. A good candidate for 8 would be the one that minimizes the variance of 
the second term in equation ( I50t . 



: arg min Var 

5>o 



Fg{S{T))W{T) 



and deriving w.r.t. 8: 



Var 



Hg{S{T))W{T) 



Var 



xaT 



W{T) {S{T)-K)-8 



xaT 

At this point we find 5 such that: 

W(T) {S(T)-K)-8 



xaT 



28 



xaT 



then 



28 

{S{T]-K)W{r 



0, 



8 = 



W{T) 
xaT 



(51) 



(52) 



(53) 



(54) 



In order to have an operative parameter we then consider the following approximation: 



Var 


'iS{T)-K)W{T)] 


[ xaT \ 


Var 







(55) 



As already mentioned, the considerations here above are still valid for the computation of the 
greeks of the options we are considering. As already illustrated, it suffices to replace 
with the Skorohod integral and that is the reason why we have always shown the term this 
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Table 1 Inputs Pai'ameters 



S,(0) 



100, Vi=l...,M 

5% 

1 



r 



T 



Pil 



10%+i^40% i=l...,M 



= 50% iJ = l...,M 



term explicitly in the calculations above. The same substitutions must be made to calculate 
each 5 for the each component of the Delta of the call type Asian basket and exotic options 
since these results hold true in the multidimensional setting as well. 

In the spirit of adaptive MC techniques (see for instance Jourdain (6)), the variance 
above can be easily estimated by a MC simulation and then, by fixing the same random 
draws, one runs a second MC simulation in order to estimate the greeks. 

In the case of one dimensional digital options the computation is slightly different. 
Kohatsu-Higa and Patterson fl] claim that a good candidate for 5 is: 



The above parameter can be easily estimated by an adaptive MC simulation in the multidi- 
mensional setting as illustrated for call-type options. 

We note that in our formulas the computation of the k-th delta depends only on the k- 
th component of the Skorohod integral making the localization technique easier to apply 
and the parameter 5 easy to calculate. Once more, we remark the fact that these variance 
and computational reduction considerations are not possible without using the Malliavin 
Calculus approach. 



5 Numerical Investigations 

In this section we discuss the results of the (R)QMC estimation based of the proposed ap- 
proaches. We consider M = 5 and M = 10 underlying securities and an equally-spaced time 
grid with A' = 64 time points. Hence, the effective dimension of the (R)QMC simulation is 
either 320 or 640. We estimate the multidimensional Deltas (with respect to each underlying 
asset) of each contract discussed before. The parameters chosen for the simulation are listed 
in Tabled 

We adopt RQMC simulations, based on the enhanced version illustrated in Section 14X1 
and consists of 32 replications each of 2048 random points. These random draws are ob- 
tained from a Matousek affine plus random digital shift scrambled version (see Matousek JOj) 
of the Sobol sequence satisfying Sobol's property A (see Sobol |19|). We also avoid gen- 
erating the 320 or 640-dimensional Sobol' sequence by using a Latin Supercube Sam- 
pling (LSS) method (Owen 1 14|). Briefly, this sampling mechanism is a scheme for cre- 
ating a high-dimensional sequence from sets of lower-dimensional sequences. For instance. 




(56) 





(57) 
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a 640-dimensional low discrepancy sequence can be concatenated from 13 sets of a 50- 
dimensional low discrepancy sequence by appropriately randomizing the ran order of the 
points (the last concatenation neglects the last 10 dimensions). For theoretical justification 
of the LSS method, see Owen |[T4l . 

The computation is implemented in MATLAB on a laptop with an Intel Pentium M, 
processor 1.60 GHz and 1 GB of RAM. We compute all the optimal columns for the LT 
technique in Section|4T| Such an LT construction is optimal if the integrand function is the 
payoff of the option and hence is optimal for price estimation. In contrast, our goal is the 
computation of the Deltas and this would not seem to be the optimal choice. However, if 
we would have applied the LT for the integrand function given by the Malliavin approach 
we would have got as many LT-decomposition matrices as the number of assets (one for 
each delta). This setting would remarkably increase the CPU time making the estimation 
less convenient. The numerical experiments below justify our assumption. 



5.1 Call with fix and floating Strike 

As a first experiment, we compute the Deltas of an Asian basket option with fixed and float- 
ing strike. We compare the estimated values of the Deltas and the accuracies obtained with 
different approaches: finite differences, localization with different parameters and finally 
localization coupled with adaptive parameters. The choice of the parameters for the local- 
ization and finite difference techniques is of fundamental importance because it influences 
the variance of the estimator (see for instance L'Ecuyer [8 |). The numerical derivative is 
often calculated assuming 5 = 1% (in our case 1% of the initial price of the underlying 
securities); this may not be the optimal choice. In addition, in the multidimensional compu- 
tation (gradient estimation) one should consider different 5. Our approach based on adaptive 
techniques overtakes this problem by calculating the parameters on the fly. These parameters 
are optimal meaning that they provide the minimal variance of the estimator (in the sense 
described in Section l4~2t . Tableland Table[3]show the results with different approaches ob- 
tained for an at-the-money Asian call with fixed and floating strike and M = 10 underlying 
assets. All the estimated values are in statistical accordance but display different accuracies. 
The finite difference errors are higher than those obtained with localization (with the exemp- 
tion of 5 = 5%). In particular, when the strike is floating, this technique returns a completed 
biased Delta associated with the highest volatility. Finally, finite difference estimations re- 
quire a computational effort that is 2.43 times higher that those obtained with localization. 
The adaptive localization and standard localization perform equally well with the former 
having slightly better precision and the advantage of selecting better localization parameters 
for each component. 

In order to have a complete picture of the sensitivity of the discussed techniques, we 
repeat the experiment considering only M = 4 assets and several strike prices. This further 
analysis cannot be performed for Asian option with floating strike. Figure |2] and [T] show the 
estimated Deltas and errors, respectively. Since for at-the-money options the finite difference 
approach provided lower accuracy, we avoided to report its results. In term of precision, in 
this setting as well, the standard localization with 5 = 1% and the adaptive localization 
return the most accurate results. In particular, these two approaches perform equally well 
with the former one having a more constant trend across all the moneyness. 
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0.41 
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0.58 
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1.40 


6.16 
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Table 2 Call Option with Fixed Strike, M = 10: At-the-Money Deltas and Errors (x 100). 
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Table 3 Call Option with Floating Strike, M = 10: Deltas and Errors (x 100). 



5.2 Digital Call 

The aim of this subsection is to describe the results of our numerical investigation assuming 
Asian digital options. The following discussion and description have a double purpose. Since 
the payoff of digital option can be seen as the derivative (in the sense of distribution) of 
the payoff of a call option, the methodology and the localization parameters described in 
Section [STI can be be rearranged and used to compute the Gamma (and cross sensitivities 
in the multidimensional setting) of a call option (naturally with some changes). In addition, 
the Delta of a digital option is a more demanding task due to the irregular payoff that is 
pathologically not differentiable. 

We repeat the organization of our discussion as done for the Asian call options and con- 
sider only a fixed strike price. Table |4] shows the estimated multidimensional Deltas and 
their errors for an at-the-money digital option on M = 10 underlying securities. The best ac- 
curacy with the standard localization technique is not achieved anymore with 5=1%, that 
means that in some situations it is not the optimal choice. In contrast, the adaptive local- 
ization is the best performing technique in terms of precision. It returns better localization 
parameters that provide an unbiased estimator with lower variance. 

As done before, we run a QMC simulation considering only M = 4 assets and analyze 
the results by varying the strike price. Figure |3]and|4]show the estimated Deltas and errors, 
respectively. Once more the adaptive localization approach displays the lowest error. 
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Fig. 1 Call Option with Fixed Strike, M = 4: Estimation Errors. 

Adaptive: Solid Line, Loc. S = 0.01: Dashed line, Loc. 5 = 0.1: Dotted Une, Loc. 5 = 0.05: Dash-dotted 
line. 
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Table 4 Digital Option with Fixed Strike, M = 10: At-the-Money Deltas and Errors. 
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Fig. 2 Call Option with Fixed Strike, M = 4: Estimated Deltas with the Adaptive Localization. 



5.3 Exotic Option 



As a last experiment we perform a QMC numerical simulation in order to estimate the Deltas 
of an exotic option. Table|5]and Figures|5]and|6]present the results of this experiment. In this 
last example all the approaches perform equally well, and the exotic structure of the payoff 
makes its estimator unsensitive to the different localization parameter. The finite difference 
is also performing well but is less precise if we take into account the computational burden 
that is 2.61 times higher. 



6 Concluding Remarks 



In this paper we have investigated the use of Malliavin calculus in order to calculate the 
Greeks of multiasset complex path-dependent options by QMC simulation. As a first re- 
sult we have derived the multidimensional version of the formulas obtained by Montero 
and Kohatsu-Higa lllll in the single asset case. The multidimensional setting shows the ad- 
vantage of the Malliavin Calculus approach over alternative techniques that have been previ- 
ously proposed. These different techniques are hard to implement and in particular, are com- 
putationally time consuming when considering multiasset derivative securities. In addition, 
their estimators potentially display a high variance (see for instance Chen and Glasserman 
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Fig. 3 Digital Option with Fixed Strike, M = 4: Estimation Errors. 

Adaptive: Solid Line, Loc. 5 = 0.01: Dashed line, Loc. 5 = 0.1: Dotted Une, Loc. 5 = 0.05: Dash-dotted 
line. 
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Table 5 Exotic, M = 10: At-the-Money Deltas and Errors (x 100). 
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Fig. 4 Digital Option with Fixed Strike, M = 4: Estimated Deltas with the Adaptive Localization. 



Ul). In contrast, the use of the generahzed integral by part formula of Malliavin Calculus 
gives enough flexibility in order to find unbiased estimators with low variance. In the multi- 
dimensional context, we have found convenient formulas that are easy and flexible to employ 
and permit to improve the localization technique. Finally, we have performed a detailed 
analysis on how the localization parameters can influence the precision of the estimators. 
Moreover, we have proposed an alternate approach, based on adaptive (Q)MC techniques 
that returns convenient parameters that can be obtained on the flight in the simulation. This 
approach provides a better precision with the same computational burden. However further 
studies would be necessary to enhance its accuracy assuming different dynamics and payoff 
functions. 

The proposed procedures, coupled with the enhanced version of Quasi-Monte Carlo 
simulations as illustrated in Sabino 1 16], are discussed based on the numerical estimation of 
the Deltas of call, digital Asian-style and Exotic basket options with a fixed and a floating 
strike price in a multidimensional Black-Scholes market. 
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Fig. 5 Exotic Option, M = 4: Estimation Errors. 

Adaptive: Solid Line, Loc. S = 0.01: Dashed line, Loc. S = 0.1: Dotted line, Loc. S = 0.05: Dash-dotted 
line. 
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